O’Hara on GitHub


knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
                      echo = TRUE, message = FALSE, warning = FALSE)

library(raster)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')

source(here('common_fxns.R'))

dir_str_mol <- '/home/shares/ohi/git-annex/impact_acceleration/stressors'
dir_str_10km  <- file.path(dir_bd_anx, 'layers/stressors_10km')

1 Summary

Collect stressor layers from CHI 2019 and reproject/aggregate them to the native resolution of the species range maps. This is a simple aggregation 11x as is done for the ocean, EEZ, MEOW, etc rasters.

1.1 Loop over all stressor layers, aggregate to spp CRS

Read in each raster, aggregate to a raster at the species projection and resolution. Because the projection is Mollweide for the source and target, and the base raster for species maps was created the same way, we can simply aggregate the original Mollweide CHI maps up by the same factor. For spp maps, the aggregation factor was 11\(\times\) to approximate 100 km2 cells: \(.934 \times 11 = 10.28 \text{ km}; \; (.934 \times 11)^2 = 105.6633 \text{ km}^2\).

Notes:

  • here we are aggregating using mean value.
  • land-based stressors will be masked to the coastal 3nmi buffer before aggregating. Real zeros in coastal zones will still be counted, but zeros farther offshore will be ignored.
  • for sea surface temperature, we will mask out sea ice prior to aggregating. The original file for the sea ice mask comes from: mazu:git-annex/Global/NCEAS-Pressures-Summaries_frazier2013/ice_mask_resampled (both .gri and .grd files needed)
### List directories that end in "final"
str_dirs <- list.files(dir_str_mol, recursive = TRUE, full.names = TRUE) %>%
  dirname() %>%
  unique() %>%
  .[str_detect(., 'final') & !str_detect(., 'archive')]
#  [1] "archive/art_fish/final"                 "archive/comm_fish/final/dem_dest"      
#  [3] "archive/comm_fish/final/dem_nondest_hb" "archive/comm_fish/final/dem_nondest_lb"
#  [5] "archive/comm_fish/final/pel_hb"         "archive/comm_fish/final/pel_lb"        
#  [7] "archive/sst_old/final"                  "art_fish/final"                        
#  [9] "benthic_structures/final"               "comm_fish/final/dem_dest"              
# [11] "comm_fish/final/dem_nondest_hb"         "comm_fish/final/dem_nondest_lb"        
# [13] "comm_fish/final/pel_hb"                 "comm_fish/final/pel_lb"                
# [15] "direct_human/final"                     "land_based/final/nutrient"             
# [17] "land_based/final/organic"               "light_pollution/final"                 
# [19] "oa/final"                               "shipping/final"                        
# [21] "slr/final"                              "sst/final"                             
# [23] "uv/final" 
### identify land-based stressors, and load coastal mask and sea-ice raster.
land_based_str <- get_str_cats() %>%
  filter(category == 'land-based') %>%
  .$stressor
coastal_str <- c('slr', land_based_str)
coastal_mask <- raster(file.path(dir_bd_anx, 'spatial/three_nm_offshore_mol.tif'))
seaice_mask <- raster(file.path(dir_bd_anx, 'spatial/sea_ice_mask_1km.tif'))

for(dir in str_dirs) {
  ### dir <- str_dirs[14] ### slr
  ### dir <- str_dirs[13] ### shipping
  message('Processing stressors in ', dir)
  files <- list.files(dir, pattern = '.tif$', full.names = TRUE)
  
  # for(f in files) {
  tmp <- parallel::mclapply(files, mc.cores = 12,
    FUN = function(f) {
      ### f <- files[1]
      str <- str_replace_all(basename(f), '_[0-9].+', '')
      str_10km_file <- file.path(dir_str_10km, 
                               str_replace(basename(f), 'mol', 'mol10km'))
      
      if(!file.exists(str_10km_file)) {
        message('Processing ', basename(f))
        str_rast_orig <- raster(f)
        
        ### check for coastal or sst stressors
        if(str %in% coastal_str) {
          message('Masking ', basename(f), ' to coastal only')
          str_rast_orig <- mask(str_rast_orig, coastal_mask)
        }
        if(str == 'sst') {
          message('Masking ', basename(f), ' to remove sea ice')
          str_rast_orig <- mask(str_rast_orig, seaice_mask)
        }
        
        str_rast_10km <- raster::aggregate(str_rast_orig, 
                                           fact = 11, fun = mean,
                                           filename = str_10km_file,
                                           overwrite = TRUE)
        ### check to make sure results match spp rasts
        # ocean_area_rast <- raster(here('_spatial/ocean_area_mol.tif'))
        # raster::compareRaster(str_rast_10km, ocean_area_rast,
        #                       extent = TRUE, rowcol = TRUE, crs = TRUE,
        #                       res = TRUE, orig = TRUE)
        ### TRUE!
    
      } else {
        # message('File exists: ', str_10km_file, '... skipping!')
      }
    })
}

1.2 Map all stressors in gifs

library(animation)

make_gifs <- function(map_stack, filename, layer_names = NULL) {
  if(is.null(layer_names)) layer_names = names(map_stack)

  x <- values(map_stack)[values(map_stack) > 0] %>%
    min(na.rm = TRUE)

  power_adj <- -1 * round(log10(mean(values(map_stack), na.rm = TRUE) * 2)) + 1
  message('Adjusting by power of ', 1/power_adj)
  map_stack_adj <- map_stack^(1 / power_adj)
  
  ramp_breaks <- c(0, seq(x, 1, length.out = 100))
  ramp_labels <- seq(0, 1, .2)
  ramp_colors <- c('grey40', hcl.colors(n = 100, palette = 'viridis'))

  capture.output({
    saveGIF({
      for(i in 1:nlayers(map_stack_adj)){
        plot(map_stack_adj[[i]], 
             col = ramp_colors,
             breaks = ramp_breaks,
             axes = FALSE,
             axis.args = list(at = ramp_labels^(1 / power_adj), labels = ramp_labels),
             main = layer_names[i])
      }}, 
      interval = 0.5, movie.name = filename, 
      ani.width = 700, ani.height = 420)
  })
  
  return(invisible(NULL))
}
str_10km_files <- list.files(file.path(dir_bd_anx, 'layers/stressors_10km'),
                             full.names = TRUE)
years <- 2003:2013

str_map_df <- data.frame(str_file = str_10km_files) %>%
  mutate(str = str_replace(basename(str_file), '_[0-9].+', ''),
         year = str_extract(basename(str_file), '[0-9]{4}') %>%
           as.integer()) %>%
  group_by(str) %>%
  filter(min(year) <= 2003 & max(year) >= 2013) %>%
  ungroup() %>%
  filter(year %in% years)

stressors <- str_map_df$str %>% unique()

for(stressor in stressors) { ### stressor <- stressors[14]
  ### Animate the results
  gif_file <- here('_setup', sprintf('figs/str_movie_%s.gif', stressor))
  
  str_maps <- str_map_df %>%
    filter(str == stressor) %>%
    .$str_file
    
  if(!file.exists(gif_file)) {
    message('creating ', gif_file)
    rast_files <- str_maps
    map_stack <- stack(rast_files) 
    make_gifs(map_stack, 
              filename = gif_file,
              layer_names = paste(stressor, years))
  }
  
  cat(sprintf('![](%s)', gif_file))
}